Phonon-assisted relaxation and tunneling in self-assembled quantum dot molecules 
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We study theoretically phonon-assisted relaxation processes in a system consisting of one or two 
electrons confined in two vertically stacked self-assembled quantum dots. The calculation is based on 
a k p approximation for single particle wave functions in a strained self-assembled structure. From 
these, two-particle states are calculated by including the Coulomb interaction and the transition 
rates between the lowest energy eigenstates are derived. We take into account phonon couplings via 
deformation potential and piezoelectric interaction and show that they both can play a dominant 
role in different parameter regimes. Within the Fermi golden rule approximation, we calculate the 
relaxation rates between the lowest energy eigenstates which lead to thermalization on a picosecond 
time scale in a narrow range of dot sizes. 



PACS numbers: 73.21. La, 73.63.Kv, 63.20.kd 

I. INTRODUCTION 

Structures composed of two closely spaced quantum 
dots (QDs) attract much attention motivated by their 
rich physical properties as well as by possible applica- 
tions in nanoelectronics or quantum computing. A ma- 
jor factor that determines the properties of such a system 
is the electronic coupling between the dots. For closely 
spaced dots, the system spectrum can be strongly af- 
fected by tunnel coupling^^^^. Optical spectra of closely 
spaced structures indeed show clear manifestations of 
such tunneling-related effects^"—. Due to strong delocal- 
ization of carrier states over the double dot structure, 
analogous to a chemical covalent bond, such structures 
are often referred to as quantum dot molecules (QDMs) 
or artificial molecules. 

The properties of such artificial molecules are also af- 
fected by phonon-related processes which are inevitable 
in a crystal environment. Such effects will limit the feasi- 
bility of building QDM-based quantum-coherent devices 
by providing a dephasing channel for both chargei2r2& 
and spinil states. Depending on the form and local- 
ization character of the wave functions, such phonon- 
assisted transitions may either take place between two 
delocalized states or involve charge redistribution when 
an electron dissipatively tunnels to a different dot. In the 
latter case, the electron spin can be conserved^, which 
can be used to control the spin state of a magnetic impu- 
rity in one of the QDs^^. Dissipative tunneling is also in- 
teresting in a two-electron configuration, where a transi- 
tion to a doubly occupied state is only possible in a singlet 
configuration. This discrimination leads, on one hand, to 
pure dephasing of singlet-triplet superpositions^'' but, on 
the other hand, might be used to speed up the proposed 
singlet-triplet measurement protocols22.. 

From the experimental point of view, dissipative car- 
rier transfer in self-assembled structures has been stud- 
ied with optical spectroscopy methods (time-integrated 



and time-resolved photoluminescence, and photolumines- 
cence excitation experiments) both in lateral double-dot 
systemsiSii^i^ as well as in stacked QpMsM'M. and QD 
chains (both stacked and lateral) '^^i'^^ . Various mech- 
anisms have been invoked to account for the observed 
properties. In most cases, the kinetics is attributed to 
tunneling^ir— In some other experiments^"—, sig- 
natures of radiative (Forster-like) transfer have been ob- 
served. Coulomb scattering^^ and thermally activated 
processes^L^i also seem to play an important role, at 
least in some systems. 

The variety of investigated structures and probable 
transfer mechanisms is reflected in a relatively wide dis- 
tribution of measured transfer times. While the transfer 
in general takes place on time scales shorter or compa- 
rable with the exciton life time (which is necessary for 
the process to be observable in an optical experiment), 
the observed times range from tens of picoseconds22i^ to 
several nanosecondaS,. In most cases, however, trans- 
fer times between hundreds of picoseconds and a few 
nanoseconds are observed. This experimental situation 
indicates that carrier kinetics in coupled QD structures 
is a rich and complex problem which most likely can- 
not be solved by proposing a unique, universal theory. 
Therefore, it seems reasonable to undertake a systematic 
theoretical study of various carrier transfer processes and 
to identify conditions in which one or another mechanism 
is expected to dominate the system properties. Such a 
theoretical analysis of individual transfer mechanisms has 
in fact already started with several works devoted to elec- 
tron tunneling (in simplified confinement models)^"— 
and some studies of the Forster-like transfer—"—. 

In this paper, we develop a theoretical description for 
phonon-assisted relaxation and charge transfer (tunnel- 
ing) in a structure composed of two vertically stacked 
quantum dots formed in the Stransky-Krastanov self- 
assembly process by strain-induced spontaneous QD nu- 
cleation in the second layer on top of the QDs formed 
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FIG. 1: (Color online) The geometry of the QDM structure. 



in the first laye r A reliable calculation of tunneling 
rates requires reasonably precise knowledge of the elec- 
tron wave functions. For a strained self-assembled struc- 
ture presently under consideration, this implies the need 
to calculate the strain and then to find the single particle 
wave functions, e.g., by a fc-p method. Then, Coulomb 
interactions can be included for a two-electron system 
within the standard configuration-interaction approach. 
The h p method for strained semiconductor heterostruc- 
turcs is a well established procedure that has been used 
for QDs, QDMs and other nanostructures^"— . Recently, 
this method has been combined with the standard ap- 
proach to carrier-phonon coupling in a study of confined 
polarons*". Here, we apply a simplified version of this 
method^* assuming a cylindrical symmetry of the struc- 
ture. This is motivated not only by economy of compu- 
tations but, more importantly, by the need to derive the 
wave functions in a form suitable for efficient calculations 
of carrier-phonon couplings and the following modeling 
of phonon-assisted relaxation. 

The paper is organized as follows. In Sec. [ITl we de- 
fine the model of the system under study. In Sec. IIIIl 
we discuss the strain fields in the structure. The one- 
electron and two-electron states in the QDM are found 
in Sec. [IV] and Sec. |Vl respectively. In Sec. IVH phonon- 
assisted relaxation for one- and two-electron states is dis- 
cussed. Concluding remarks and discussion are contained 
in Sec. IVlll 



II. MODEL 

We consider a QDM formed by two self-assembled InAs 
dots in a GaAs matrix. The geometry of the structure 
as used in our modeling is shown in Fig. [TJ The QDs 
are modeled as two spherical segments with base radii 
ri, r-i and heights ll2, respectively. Throughout the 
paper, the aspect ratio of the two dots will be held con- 
stant, Hx/n = H2/r2 = 0.37. Both dots are placed 
on a wetting layer with thickness 7?wl- The dots are 
separated by a distance D (base to base). A diffusion 
layer of a very small thickness -ffdiff = 0.3 nm is included 



at the contact between the two materials, in which the 
InAs concentration varies linearly. Apart from this, the 
InAs content is assumed to be 100% inside the dots and 
the wetting layers and outside. The parameters of the 
modeled structure are collected in Table [D 

Our model includes the case of a single electron in 
the QDM as well as of two electrons coupled by the 
Coulomb interaction. The carriers interact with bulk 
acoustic phonons via standard deformation potential and 
piezoelectric interaction mechanisms. 

The modeling proceeds in three steps: (1) Determina- 
tion of the strain distribution; (2) Calculation of the wave 
functions for single- and two-electron states; (3) Calcula- 
tion of relaxation rates. As each of these steps involves a 
specific formalism, the corresponding details of the model 
will be subsequently introduced in the following sections. 



III. STRAIN 

In this section, we calculate the strain present in the 
inhomogeneous structure. 

The strain fields in the system will be described by the 
strain tensor 



dui{r) duj{r) 



where u{r) is the displacement field at the point r in the 
crystal. The elastic energy of the inhomogeneous system 
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Here Cij are position-dependent elastic constants (see 
Table H] for the values), a = (Cn + 2Ci2)(ai/aG - 1) 
in a QD and a = in GaAs, where ai and are the lat- 
tice constants of InAs and GaAs, respectively. The last 
term in Eq. ^ accounts for the mismatch of the lattice 
constants, shifting the equilibrium of the InAs crystal 
lattice to the state appropriately stretched with respect 
to the ideal InAs crystal. Since the strain is calculated 
with respect to the GaAs lattice and GaAs crystal coor- 
dinates are used, after minimizing the strain energy the 
results for the InAs dots must be rescaled to yield phys- 
ical strain, according to^ 
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£■ ■ = ( 
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For an axially symmetric structure, it is convenient 
to perform the computation in cylindrical coordinates 
(p, (/>, z). Therefore we denote the components of the dis- 
placement in the local reference frame as Up,u^, Uz and 
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GaAs 


InAs 




Lattice constant 


a 


0.56532 nm 0.60583 nm 




Elastic constants 


Cii 


12.11 • 10^° N/m^ 8.329 ■ 10^" 


N/m^ 




Cl2 


5.48 ■ 10^° 


N/m^ 4.526 ■ 10^" 


N/m^ 




C44 


6.04 ■ 10^" 


N/m^ 3.959 ■ lO"' 


N/m^ 


Band structure parameters 


Eco 


0.95 eV 









-EvO 


-0.57 eV 


-0.42 eV 






A 


0.34 eV 


0.43 eV 






Po 


9.89 eV 


9.19 eV 




Deformation potentials 




-9.3 eV 


-6.66 eV 








0.7 eV 


0.66 eV 






b 


-2.0 eV 


-1.8 eV 




Speed of sound 


Cl 
Ct 




5150 m/s 
2800 m/s 




Crystal density 


Q 




5300 kg/m^ 




Piezoelectric constant 


dp 




-0.16 C/m^ 




Relative dielectric constant 






12.9 





TABLE L System parameters used in the calculations. 



define the corresponding components of the strain tensor, 

" dp ' 
_ 1 / du^ 
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1 / dup duz 



2\dz dp J' 
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2\dz p 



We will look for the minimum of E^i in the class of ax- 
ially symmetric displacement fields, that is, u,p — and 
dur/dcj) — duz/d4> = 0. With this assumption, the inte- 
gration over (j) in Eq. ^ can be performed analytically 
and one gets 

E,i = TT J dpp J dz[Cnel,+D{elp + el^)+4C44e% 



-2a(epp -I- e, 



0], 



(3) 



where D = 3Cii/4 -I- C12/4 -I- Cm/2 and F = Cii/4 + 
3C12/4-C44/2. 

The displacement field minimizing Ed is found by 
the conjugate gradient method on a square grid of 1000 
points along z and 666 point along p, representing a cylin- 
der with the height of 60 nm and the radius of 40 nm. The 
boundary conditions represent displacements due to the 
equilibrium strain in a system with two wetting layers, 
which corresponds to the actual situation at large dis- 
tances from the dots. A combination of discretizations 
with forward and backward representations of derivatives 



is used to avoid discretization-induced oscillations^. As 
an example of the result, a strain map representing the 
hydrostatic and axial strain across the structure for a 
selected geometry is shown in Fig. O 



IV. SINGLE ELECTRON STATES 

In this section, we calculate approximate wave func- 
tions for a single electron confined in the nanostructure. 
This is done within a variational multi-component enve- 
lope function scheme^ based on the fact that the confine- 
ment volume is large compared to the crystal lattice cell 
and that the local system parameters change relatively 
slowly on atomic scales. In this approach, one finds the 
values of effective masses and band edges at a given point 
by solving the bulk k p model with strain and composi- 
tion equal to those present at a given point. This yields 
the band edge position, which is used as the local ef- 
fective potential, as well as the band curvatures, which 
define the components of the effective mass tensor at a 
given point of the inhomogeneous heterostructure. 

The conduction band structure in a strained system 
is determined from the 8-band k p (Kane) Hamiltonian 
with strain-induced terms (Bir-Pikus Hamiltonian) using 
the Lowdin elimination^. The part of the Hamiltonian 
coupling conduction and valence band states is^ 



i?c-v = |ct) [-VSV^hhtl - U (\/2(lht| - (soti) 

-v^((ih;| - V2(so;|) 

(so II 



((lht| + V2(sot 



-^H.c. 
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FIG. 2: (Color online) The strain distribution in the structure 
for D = 12 nm, r\ = 10 nm, r2 ~ 10.5 nm, Hi — 3.7 nm, 
H2 = 3.885 nm. In (a), the hydrostatic strain is shown; in 
(b), the ajdal strain e^z — {^pp + e^^)/2 is shown. 



where 'e' 'Ih', 'hh', and 'so' denote the electron, heavy- 
hole, hght-hole and spin-orbit split-off subbands, f and I 
represent the spin orientation in a given subband. 



and 



eyj)kj 



where Pq is proportional to the interband momentum 
matrix element (see Table U for parameter values) . The 
conduction band part of the Hamiltonian is 



Hr = 



E. 



cO 



2rno 



arh 



(|et)(et| + |ei)(ei|), 



where TOq is the free electron mass, E^q is the conduction 
band edge in a bulk unstrained crystal, ac is the con- 



duction band deformation potential, and h = Tre is the 
hydrostatic strain. 

As we are interested in the corrections to the con- 
duction band energies up to and the off-diagonal ele- 
ments U and V are proportional to k we only need the 
conduction- valence band energy difference at /c = 0. In 
this limit, the diagonal terms for the valence band states 



Ehh 
Eih 

E^n 



E^o-p- q, 

EvO-p + q, 
i^vo - A - p, 



where E^q is the valence band edge of an unstrained 
crystal, A is the spin-orbit split-off parameter of an un- 
strained crystal, p = a^h, 



1 



(epp + £00) 



and av,6 are valence band deformation potentials. The 
values of material parameters are given in Table HI 

Neglecting the strain-related terms in U and V , which 
are much smaller than the purely kinetic ones, we get the 
conduction band energy up to the 2nd order in fc, 



E{k) = EcQ + Qch + 
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h'k 



2m± 



where the in-plane and z components of the effective mass 
tensor are 



2AE, 
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Ev 
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2Ev 



6AEi^ 



3AK, 



Ep 



3AEih 3AEs, 



E„ for 



where Ep = 2moP§/h'^ and AEi = Eco + ach 
i = hh, Ih, so. Note that AEi are position dependent. 

The dynamics of an electron in the strained nanostruc- 
ture in the present approach is defined by the conduction 
band edge at a given point, Ec{p,z) = Eco + ach{p,z) 
(which depends on the local strain), and on the effec- 
tive masses, which also vary across the structure. Fig. [3] 
shows an example of the profiles of the conduction band 
edge as a function of z at three different values of p. In 
Fig. 21 we show the spatial maps of the radial and axial 
components of the electron effective mass. Some strain- 
induced anisotropy of the effective mass can be seen. The 
value of the radial component, m± sa 0.05mo is close to 
the bulk GaAs value and much higher than the bulk InAs 
value of 0.023mo. It is roughly constant within the vol- 
umes of the two dots. The axial component is lower 
(about 0.04mo) and shows some gradient along the QD 
axis, with higher values towards the top. 
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FIG. 3: (Color online) The band edge profiles along z for 
three fixed values of p = 0.3 nm (red solid line), 3 nm (blue 
dashed line), and 8 nm (green dotted line) for the structure 
as in Fig. [2] 



The envelope function of an electron is found from the 
Schrodinger equation with the Hamiltonian 



H ^ - 



d 



d d 



d 



dx 2m±{p, z) dx dy 2m±{p, z) dy 
d d 



dz 2mz{p, z) dz 



Ec{p,z) 



Following the concept of 'adiabatic' separation of 
variables^i, we first numerically solve the one- 
dimensional equation along the strongest confinement di- 
rection at each p, 



d 



d 



dz 2mz{p, z) dz 



X{p,z) = E{p)x{p,z). 



The lowest two solutions to this equation, Xi,2(p, z), rep- 
resent the lowest subband of confined states in the dou- 
ble well system. The corresponding two branches of p- 
dependent eigenvalues, Ei 2{p) can be interpreted as ef- 
fective potentials for the radial problem. 

Next, we apply the Ritz variational method^, looking 
for the stationary points of the functional 

F[i;] - (^liJIV') 

in the class of normalized ansatz functions 



(4) 



where M is the angular momentum. Upon transforming 
to cylindrical coordinates and imposing the normaliza- 
tion via the Lagrange multiplier A, we write the func- 
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FIG. 4: (Color online) The components of the electron ef- 
fective mass for the structure as in Fig. (2] (a) the radial 
component, m±; (b) the axial component, m^. 



tional F[7ji\ in the explicit form 



2m±{p, z) 



E^{p) + — 

P 



Y] / Pdp^lip) 
„. Jo 

POO 

Y] / Pdpfl{p)fa{p) - 1 



'Pa{p) 



We discretize the functional on the same lattice 

that was used in the computation of the strain. As 
the functional is quadratic, the stationarity requirement 
with respect to the values at the discrete points can eas- 
ily be cast into the form of a matrix eigenvalue prob- 
lem for the components {Lpi,ip2). By virtue of the Ritz 
theorem^, the corresponding eigenvalues A„, n = 0, 1, . . . 
approximate the energy eigenvalues of the original prob- 
lem, while the eigenvectors, representing the components 
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FIG. 5: (Color online) (a-c) The single-electron energy levels 
for three structures with a fixed size of the lower dot as a 
function of the size (base radius r2) of the upper one for three 
dot separations as shown. Here ri = 10 nm, Hi = 3.7 nm, 
H2/r2 = 0.37. The energy reference level is 0.8 eV above 
the conduction band edge of unstrained bulk InAs. (d-f) The 
corresponding probabilities piow of finding the electron in the 
lower half of the system as a function of the size of the upper 
dot, in the ground state (labelled "g") and in the first excited 
state (labelled "e") of the system. 



[yyi^\p) , 'P2^\p)] at the discrete lattice points in the ra- 
dial direction, are used to construct the electron eigen- 
functions ^nii") = ipnip, z, (p) according to the ansatz for- 
mula In this paper, the discussion will be restricted 
to the two lowest states for Af = 0, corresponding to the 
tunnel-split ground state of the double-dot system. 

The single particle eigenenergies found within our ap- 
proach for three different distances D between the dots 
are shown in Fig. [5f a-c). In these calculations, the shape 
of the lower dot is kept constant, while the base radius 
r2 of the upper dot and its height H2 are varied, with 
Hijri constant. Electronic (tunnel) coupling between 
the dots leads to the appearance of an anticrossing pat- 
tern near the point where the dots become equal. The 
width of the anticrossing is 2t, where the phenomenolog- 
ical "tunnel coupling" i is very well fitted by the formula 
Ini/io = —kD, with k = 0.59 nm^^ and to ~ 1.58 eV. 
The value of k is consistent with the one-dimcnsional 
semiclassical formula k = y^2mz{V — E)/h if one uses 
ruz — 0.062mo as found in the area between the dots 
(see Fig. 2]), the potential barrier height V = 1010 meV 
(Fig. [3]) and the electron energy E = 800 meV (Fig. [5]). 

In accordance with the spectral anticrossing, the elec- 
tron occupations for the ground and first excited states 
are transferred between the dots, as shown in Figs. El^d- 
f). The exact resonance point, where the two occupations 
are equal to 1/2 for both states (corresponding to de- 
localized symmetric and antisymmetric wave functions), 
appears for r2 slightly smaller than ri which results from 
the strain field in the absence of mirror symmetry in the 



structure. 

The procedure proposed here, involving the variational 
problem for a two-component envelope, allows for mixing 
of the two manifolds of states related to the two functions 
Xi(p, z) and X2(p, z), which is essential when the two QDs 
are of similar size or when the thinner dot has a larger 
in-plane size, so that a crossing of the one-dimensional 
solutions appears at a certain value of p. 



V. COULOMB INTERACTION AND 
TWO-ELECTRON STATES 

We will find the two-electron states in the restricted 
basis of low-energy configurations of the two electron sys- 
tem. We discuss the situation when the energy differ- 
ence between the ground states in the two dots is smaller 
than the intra-dot excitation energy (the latter is about 
50 meV) . Then the two lowest single-particle states found 
in Sec. IIVI correspond to an electron in the ground state 
of one of the dots or, near the resonance, to a delocalized 
superposition of the two ground states. 

Let Qn.a, cin a dcuote the annihilation and creation op- 
erators for an electron in the state n = 0, 1 with the wave 
function ipn{f) and spin a. The low-energy two-electron 
configurations split into one triplet state (of no interest 
in the present discussion) and three singlet states 



|0) = aj^4_^|vac), 



(5a) 



ID = ^^li±ili|vac), (5b) 



V2 

1 2) = a|^aj^ I vac). 



(5c) 



where |vac) is the vacuum (empty dot) state. 

The Hamiltonian of the interacting two-electron sys- 
tem has the form 

H enai^ana + ^YY1 '"^jkiala^^lcr'^ka'aia, (6) 



n,<7 



ijkl a^a' 



where 



Vijkl 



1 



jrMr')Mr). (7) 



Here e is the electron charge, Sq is the vacuum permit- 
tivity, and El- is the dielectric constant of the semicon- 
ductor. Some technical details concerning the calcula- 
tion of Coulomb matrix elements for the wave functions 
obtained within the variational two-component envelope 
function scheme in Sec. IIVI are given in the Appendix. 

In Fig. [51 we show the three lowest spin-singlet eigen- 
states of the interacting two-electron system as a function 
of the size of the upper dot with the lower dot kept fixed. 
The central resonance occurs when the dots are close to 
identical and involves the doubly occupied configurations 
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FIG. 6: (Color online) (a-c) The two-electron energy levels 
for three structures with a fixed size of the lower dot as a 
function of the size (base radius r2) of the upper one for three 
dot separations as shown. Here r\ = 10 nm, H\ = 3.7 nm, 
Hijr^ = 0.37. The energy reference level is 1.6 eV above the 
conduction band edge of unstrained bulk InAs. In (b), the 
electron configurations corresponding to the spectral branches 
are shown, with the first and second digits corresponding to 
the number of electrons in the upper and lower dot, respec- 
tively, (d-f) The corresponding average numbers of electrons 
in the lower half of the system as a function of the size of the 
upper dot, in the ground state (labelled "g") and in the first 
and second excited states (labelled "el,e2") of the system. 



(0, 2) and (2, 0) which, at the resonance point, have sim- 
ilar energy. This anticrossing is very narrow (less than 
0.1 meV for D — 12 nm) since the two states involved 
differ by the location of both electrons [see Figs.[ni[d)-(f)] 
and, therefore, are coupled only by very small exchange- 
like Coulomb terms. Only for the smallest inter-dot dis- 
tance considered, D = 9 nm, this splitting becomes larger 
due to stronger mixing of configurations and incomplete 
electron localization in the two states (which allows the 
configurations to be coupled by single-electron tunnel- 
ing). The other two anticrossings occur at the degener- 
acy point between the singly occupied (1,1) configura- 
tion (favored by the Coulomb repulsion) and the (0, 2) 
or (2,0) configuration with two electrons in the larger 
dot. One can notice that these two splittings are wider 
than those appearing between the single-electron states, 
shown in Fig. [5] (for instance, 2 meV vs. 1.5 meV for 
D = 12 nm). This is due to the fact that the anticrossing 
of two-electron configurations is enhanced by Coulomb 
termsi^. 



(8) 

where the coupling constants Fs,nn' (?) have the symme- 
try Fs^nn'iq) = Ps,n'ni^^)- ^hc intcr-lcvel energy dis- 
tance in our structure is smaller than the optical phonon 
energy. Therefore, only acoustic phonons are relevant in 
our model. We include the deformation potential (DP) 
coupling to longitudinal acoustic (LA) phonons and the 
piezoelectric (PE) coupling to LA as well as transverse 
acoustic (TA) phonons. The coupling constant for the 
DP coupling mechanism is given by 



i^iw(g) 



/ hq 
2pVci 



(9) 



where g is the crystal density, V is the normalization 
volume of the phonon modes, ci is the longitudinal speed 
of sound (see TableUfor parameter values), and the form 
factor is defined as 



-F n 



(10) 



The coupling element for PE interactions reads 



dpe 



2pVcsq eo£i 



Ms{q)TnrAq), (11) 



where Cg is the speed of sound (s = l,t denotes the LA 
and TA phonon branch, respectively) and dp is the piezo- 
electric constant. The function Ms{q) does not depend 
on the value of the phonon wave vector, but only on its 
orientation. For a zinc-blende structure, it reads 



Ms{q) = qx [ies,q)yqz + {es,q)zqy] 
+Qy [{es,q)zqx + ies,q)xqz] 

+ qz [{es,q)xqy + ies,q)yqx] 



(12) 



where is.q is the unit polarization vector for the phonon 
wave vector q and polarization s, and q = q/q. We 
choose the following phonon polarization vectors 

e\^q = q = (sin 6* cos sin 6* sin 0, cos 6*) , (13) 
eti.q = (- sin0, cos(^, 0) , 
et2,<j = (cos cos cos ^ sin 0, — sin 0) , 

for which the functions Mgiti) read 



VI. PHONON-ASSISTED RELAXATION 

In this section, we discuss the phonon-assisted relax- 
ation between the single-electron states and between the 
two lowest two-electron states. 

The coupling between the electrons and phonons is de- 



Mi(q) = -sin(26')sin6'sin(2(/)), 

Mti(q) = sin(26')cos(2(/)), 

Mt2{q) = sin6'(3cos^6'- l)sin(2(?!)). 



(14) 



In what follows, we will assume that higher states are 
separated by an energy much larger than k^T ^ where 
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fee is the Boltzmann constant and T is the temperature. 
Then, the kinetics leading to thermahzation of the occu- 
pations of the two relevant levels can be characterized by 
the occupation of the upper state, 



-7* 



where no is the initial occupation, 7 is the relaxation 
(thermahzation) rate and 



1 



''cq 



gA_E/(feE 



1 



is the equilibrium occupation, where /S.E > is the en- 
ergy separation between the two states. 

Thus, given the initial condition and the energy dif- 
ference Ai?, the thermahzation kinetics is determined by 
the relaxation rate 7 (or the relaxation time r = 7^^) 
which will be found in the following sections, first for a 
single-electron, then for the two-electron case. 



A. Single electron relaxation 

For a single electron system, the thermahzation rate 7 
can be found directly from Eq. ([5]) using the Fermi golden 
rule. The result can be written in the form 



7 = 27r [2nB(A£;) + 1] J{AE/h), 



(15) 



where 



1 



^AE/{kBT) _ I 



is the Bose distribution and the spectral density J{Ld) is 
given by 

•^(^) = P E \FsMQ)\^^i^ ~ ^<i,s), (16) 



where Fsfii{q) is the total coupling for the branch s, 
that is, Fifli{q) = Fl^f\q) + Fj^Piq) and Fs,oi(q) = 

^tof'l?) 1*^^ ^ — tl,t2. In fact, due to different parity 
of the DP and PE couplings (as functions of q) the two 
contributions do not interfere and the spectral density 
(hence, also the thermahzation rate) can be split into 
the corresponding two parts J^^^^{uj) and j'^^^^w). 

In order to find the thermahzation rate, we calculate 
the form-factors defined in Eq. (fTO|) using the single- 
electron wave functions found for the strained double dot 
structure in Sec. IIVI (see Appendix). From these, we find 
the coupling constants given by Eq. ([9]) and Eq. (fTTj) and 
the corresponding spectral densities given by Eq. (|16p . 
The rate 7 then follows from the Fermi golden rule for- 
mula, Eq. (fT5|l . 

The single particle relaxation rates are shown in 
Figs. EJa-c) as functions of the upper dot size for three 
values of the inter-dot spacing (for the same sample ge- 
ometries as in Fig. [5]) and at three different temperatures. 




(e) / 


y D = 12 nm 







10 
To (nm) 



10.3 




10.2 



(nm) 



To (nm) 



FIG. 7: (Color online) (a-c) Thermalization rate for one elec- 
tron states in a structure witli ri = 10 nm, Hi = 3.7 nm 
and H2/r2 ~ 0.37 for tiiree different inter-dot separations at 
T = OK (red solid line), 20 K (blue dashed line), and 40 K 
(green dotted line) . (d-f ) Contributions to tiie thermalization 
rate from the DP coupling (blue dashed line) and PE coupling 
(green dotted line) as well as the total rate (red solid line) at 
T = K. 



These three plots show that both the magnitude and the 
size dependence of the relaxation rate is different in these 
three cases. The interpretation of this behavior can be 
based on the Fermi golden rule in the form of Eq. (|15p . 
where the essential role is played by the spectral density 
defined in Eq. (HH) and plotted (for D = 12 nm) in FigM 

The overall magnitude of the spectral density de- 
pends on the spatial overlap between the wave func- 
tions corresponding to the states involved in the transi- 
tion. It is, therefore, large at the resonance and becomes 
smaller as the system is shifted off the resonance point 
Fig-ISUa). Apart from this, the spectral density shows os- 
cillations in its high-frequency tail which are due to the 
essentially one-dimensional emission of short wavelength 
phonons along the strongest confinement direction^^ . As 
we deal with two confinement centers displaced along 
the same direction, interference effects appear and the 
phonon emission amplitude has a maximum whenever 
oj = (2n -|- 1)ttc/ D for an integer n. Moreover, the en- 
velope of the spectral density decays at high frequencies 
since the short wave length phonons are not effectively 
coupled to the relatively weakly confined electron states. 

In the case of closely stacked dots [Fig.[7Ka)], the tun- 
nel splitting of the QDM electron states is large and the 
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fi(B(meV) /i»(meV) 



FIG. 8: (Color online) (a) The spectral density of the phonon 
reservoir as a function of frequency for three different values of 
the upper dot radius in the resonance area, red solid line: r2 = 
9.96 nm (exact resonance); blue dashed line: r2 = 9.9 nm; 
green dotted line: r2 ~ 97.8 nm. (b) The contributions to 
the spectral density at the resonance from the deformation 
potential coupling to LA phonons and from the piezoelectric 
coupling to LA and TA phonons. 



frequency of the emitted phonons always lies far in the 
tail of the spectral density. This is reflected by the very 
low relaxation rate. The oscillations of the spectral den- 
sity are clearly marked in the values of the relaxation 
rate. When the dots are separated by a larger distance 
[Fig.[7^b)], the resonance becomes narrower and now the 
resonant frequency lies in the region of large spectral den- 
sity. When moving away from the resonance, the relax- 
ation rate drops down primarily due to the decreasing 
overlap of the wave functions. This leads to a narrow 
peak in the dependence of the relaxation rate around 
the resonance value. Still, oscillations are visible in the 
slopes of this peak. For even larger inter-dot distances 
[Fig-IZIc)], the resonance becomes very narrow. Corre- 
spondingly, the overlap between the wave functions de- 
cays almost completely already when the size of the upper 
dot is changed by a fraction of a nanometer from the res- 
onant value. Therefore, the relaxation rate is large only 
in a very narrow region around the resonance. The rates 
are also generally lower than in the previous case, which 
results from the dependence of the spectral density at 
low frequencies (~ for the DP coupling and ^ lu^ for 
the PE coupHng). 

The interplay between the shape and magnitude of 
the spectral densities for different coupling mechanisms 
[Fig. Eljb)] and the electron energies near the resonance 
is reflected also by the different contributions from the 
DP and PE couplings to the total relaxation rates. As 
can be seen in Figs. Eljd-f), the DP coupling dominates 
for large energy splittings. The reason is that this cou- 
pling is isotropic and involves LA phonons which have 
higher energies. On the contrary, the piezoelectric cou- 
pling is anisotropic and, according to Eq. p^. is sup- 
pressed for emission along the z direction that is pre- 
ferred at high frequencies. The situation changes at low 
energy splittings where the low-frequency properties of 
the spectral density are relevant. As the spectral density 
for the piezoelectric coupling decreases at w — )■ more 




D (nm) D (nm) 

FIG. 9: (Color online) (a,b) Energy splitting between the 
two lowest single-electron states as a function of the distance 
D between the dots for ri = 10 nm and r2 as shown. (c,d) 
Thermalization rate for one electron states in a structure with 
ri = 10 nm. Hi = 3.7 nm and H2/r2 = 0.37 as a function of 
D for the two values of r2 at T = K (red solid line), 20 K 
(blue dashed line), and 40 K (green dotted line). 



slowly than that corresponding to the DP coupling the 
PE coupling is the dominating mechanism in the case of 
narrow anticrossing, as can be seen in Fig.[7^f). For very 
low frequencies, all the contributions to the spectral den- 
sity are small, hence the phonon-assisted relaxation pro- 
cess becomes ineffective for small energy splitting. This 
is manifested by a dip in the thermalization rate at the 
exact resonance for D = 12 nm [Fig. [7]Jf)]. 

In Fig. [5] we show the energy splitting between the two 
lowest energy levels and the corresponding values of the 
thermalization rates 7 = t^^ as a function of the inter- 
dot distance D for two system geometries: slightly differ- 
ent dots [Fig.[i;a,c)] and identical dots [Fig.[i;b,d)]. The 
values of the rates show oscillations, resulting from the 
variation of the energy level splitting and corresponding 
to the oscillations of the spectral density, as discussed 
above. The maximum value is quite large and corre- 
sponds to relaxation times of several picoseconds, which 
results from the relative proximity of the resonance (iden- 
tical dots) in both presented cases. The maximum goes 
down and shifts to lower distances as the dots become 
different. At large distances the relaxation becomes in- 
efficient in any case. In an attempt (not shown) to com- 
pare the decrease of the rates at large D with an expo- 
nential law (as observed, at least approximately, in some 
experiments^ ^^^'^^ ), we have found a roughly exponen- 
tial decay with a coefficient consistent with the value of k 
found in Sec. IIVI This decay is, however, strongly mod- 
ulated by oscillations. This results from a small energy 
scales in our model which is comparable to strain-related 
effects as the dots are moved with respect to each other. 
This is visible in Figs.[ni[a,c), where the energy level sep- 
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aration does not tend to a constant asymptotic value at 
large D as would be expected for a simple model of two 
potential wells with fixed shapes. 



B. Relaxation in two-electron systems 

In this section, we calculate the transition rates for 
phonon-assisted relaxation between two-electron states 
l^'i), obtained from the diagonalization of the Hamil- 
tonian ([5]) in the restricted basis formed by the states 
|0), |1) and 1 2) [Eqs. ((5a| -fc)]. We first project the carrier- 
phonon Hamiltonian ^ onto the two-electron subspace, 

HfX^ = |*.)(%| G,.,(g) (fe,,, + bl 
where the coupling constants 

nm,cr 

are found based on the numerical results for the states 
jvPi). We restrict the discussion to transitions between 
the two lowest states |\I/o) and l^*!), separated by an en- 
ergy splitting AE. In the Fermi golden rule approxima- 
tion, the rate for the relaxation of the two occupations 
to equilibrium is 

7(2) = 27r [2nB{AE) + 1] J^^\AE/h), 

where the spectral density J*-^' (uj) is given by 

q,s 

The inverse relaxation times = 7^^^ resulting from 
these calculations are presented in Fig. [TOl Like in the 
single electron case, the energy level splitting in the case 
of relatively closely spaced dots {D = 9 nm) is very large 
and the resonance is very broad which results in very 
long relaxation times which do not vary considerably over 
the parameter range studied [Fig. [TUf a)]. For such high 
transition energies, only LA phonons contribute to the 
process via DP couphng [Fig. [TUTd)]. At larger inter- 
dot distances, the transition rates become large around 
the resonance points corresponding to the anticrossing of 
(1, 1) and (2, 0) or (0, 2) configurations. The structure of 
the relaxation rate as a function of the upper dot diam- 
eter r2 is similar to that discussed in the single electron 
case above. Also the relative contributions form differ- 
ent coupling mechanisms behave in the same way, with 
the piezoelectric coupling dominating at low energies. In 
general, the relaxation rates are similar to those found 
in the single electron case, since both these processes 
are physically very similar. In both cases, the electron 
tunnels between the dots and simultaneously emits one 
phonon. The only difference is that in the single electron 
case it tunnels towards an empty QD, whereas in the two 
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FIG. 10: (Color online) (a-c) Thermalization rate between 
the two lowest two-electron states for three different inter-dot 
separations at T = K (red solid line), 20 K (blue dashed 
line) , and 40 K (green dotted line) . (d-f ) Contributions to the 
thermalization rate from the DP coupling (blue dashed line) 
and PE coupling (green dotted line) as well as the total rate 
(red soUd hue) at T = K. 



electron case, there is already another electron. This ba- 
sically leads to shifts (due to Coulomb interaction) of the 
parameter regimes where the relaxation is most efficient 
from the region of identical dots to the asymmetric situ- 
ation where the difference of confinement energies com- 
pensates for the Coulomb repulsion. A similar conclusion 
has been reached in the case of gated QDM structures 
modeled by Gaussian potential wellsi^. 



VII. CONCLUSION 

We have studied phonon-assisted relaxation (thermal- 
ization) for single-electron and two-electron configura- 
tions in self- assembled quantum dots. In order to de- 
scribe the electron states in a strained structure in a 
possibly realistic (but still relatively simple) way and to 
reliably model the effect of the system geometry we have 
developed a generalized, multi-component envelope func- 
tion formalism based on the variational principle. 

Our results show that the single phonon relaxation is 
very efficient in an extremely narrow range of relative QD 
sizes near the anticrossings of energy levels but only for 
systems with a sufficiently large inter-dot distance (sev- 
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eral nm). In this case, the relaxation times can be as 
low as 1 ps, both in the single-electron and two-electron 
cases. The range of efficient relaxation becomes narrower 
as the dots are more distant from each other. Both cou- 
pling channels, piezoelectric and deformation potential, 
are important for the overall relaxation rate. The former 
dominates at low (sub-meV) transition energies. 

When the distance between the dots becomes smaller 
than about 10 nm, the energy level splitting becomes too 
large to be spanned by a single acoustic phonon (but 
still to small for an optical one). In this range of closely 
stacked dots, the tunneling times increase by orders of 
magnitudes and take values in the nanosecond range. 
For such small inter-dot distances, the energy splitting 
between the two lowest states is dominated by tunnel 
coupling and depends weakly on the size difference. As 
a result, the efficiency of the relaxation process remains 
nearly constant over a wide range of dot sizes. One can 
expect, however, that two-phonon processes^ can be im- 
portant in this range of parameters, in particular for en- 
ergy splittings exceeding the optical phonon energy. In 
general, decoherence in such systems may be dominated 
by pure dephasing due to occupation-conserving phonon 
scattering^^. 

Our findings seem to be consistent with the general 
features of experimental observations. The size range 
where the relaxation is very efficient (on picosecond time 
scales) is extremely narrow and does not exceed a few 
Angstrom, which is comparable to the lattice constant 
of GaAs. This means that such an efficient relaxation 
between the two lowest states in self-assembled quantum 
dot molecules is a rather rare phenomenon which occurs 
only for very finely tuned (accidentally or intentionally) 
dots and is unlikely to be observed in a typical sample. 
Therefore, we conclude that relaxation times on the or- 
der of at least hundreds of picoseconds should be typi- 
cal. The coupling between the dots decreases exponen- 
tially with the distance between them which reduces the 
overlap between the wave functions. Therefore, phonon- 
assisted tunneling for a spontaneously formed pair of non- 
identical dots should become inefficient as the spatial sep- 
aration between between the dots grows beyond a certain 
distance, as is indeed observed in experiments^ '^^'^^ . In 
the case studied here, that is, single-phonon relaxation 
between states separated by a few meV in energy, the 
relaxation rates undergo oscillations as functions of the 
geometrical parameters due to a structured nature of the 
phonon reservoir and the resulting interference effects. 
One should note, however, that most of the available ex- 
perimental data correspond to systems which much larger 
energy splittings. 

A more quantitative comparison is possible in the case 
of the measurements presented in Rcf. 28. Here, electron 
tunneling (that is, a transition between spatially direct 
and indirect exciton states) has been studied for a QDM 
with a fixed 10 nm spacing and energy level difference of a 
few meV, which corresponds more closely to the physical 
situation of our model. Our calculations for such param- 



eter range yield transfer times in the range of hundreds of 
picoseconds, which reasonably agrees with the measured 
time of 0.5 ns (note that a slightly different material sys- 
tem was used in that experiment and that some details of 
the system geometry are not exactly known). It will be 
interested to include the electric field in our model and to 
seek a closer correspondence with the experiment, which 
is planned as a future work. 
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Appendix A: Form factors and Coulomb matrix 
elements 



In this appendix, we briefly summarize the method 
of calculating the form factors and Coulomb matrix ele- 
ments based on the wave functions obtained within the 
variational multi-component envelope function formalism 
using the simplification offered by a cylindrically sym- 
metric system. 

Using the identity 



1 



1 



r — r' 



(2^) 



6 I qz 



one writes the Coulomb matrix element Vijki given by 
Eq. ^ in the form 



Vijkl 



(27r)3£oe 



where the form factors are given by Eq. (|TC 

We will use cylindrical coordinates for the vector r = 
[pcos (/)', /9sin0', z\ and spherical coordinates for the vec- 
tor q = [gj^ cos (/), g_L sin 0, q^], where q±^ = qsm9 and 
Qz ~ qcos9. For wave fmictions in the form given in 
Eq. ([U, one has 



J^kjiq) = Tkjiq±,qz)i 



Mj-Mki{Mj~Mk) 



(A2) 



where 



a,0 ' 

xJm,-mA<1±p)- 



(A3) 



In Eqs. (jA2[) and (jA3p . we denoted the angular momenta 
of the two states by Affc, Mj, used the identity 



1 



dd)' gilim' -m}4>' +a cos 
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where Jra is the m-th Bessel function, and introduced Substituting Eqs. (|A2[) and (jA3p into Eq. (jAip and 
the quantities integrating over </> one finds 



oo 



J-oo . '(27r)2eoe7o 

which are calculated by fast Fourier transform on the x J"ii(gsin^, gcos6')J"j.j ((7sin0, gcos^). 

grid. 
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